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We simulate numerically the full dynamics of Faraday waves in three dimensions for two 
incompressible and immiscible viscous fluids. The Navier-Stokes equations are solved 
using a finite-difference projection method coupled with a front-tracking method for the 
interface between the two fluids. The critical accelerations and wavenumbers, as well as 
the temporal behaviour at onset are compared with the results of the linear Floquet 
analysis of Kumar & Tuckerman (J. Fluid Mech., vol. 279, 1994, p. 49). The finite- 
amplitude results are compared with the experiments of Kityk et al. (Phys. Rev. E, vol. 
72, 2005, p. 036209). In particular, we reproduce the detailed spatio-temporal spectrum 
of both square and hexagonal patterns within experimental uncertainty. We present the 
first calculations of a three-dimensional velocity field arising from the Faraday instability 
for a hexagonal pattern as it varies over its oscillation period. 



1. Historical introduction 

The Faraday experiment consists of shaking vertically a container holding two immis- 
cible fluids (the lighter of which can be air) thereby inducing oscillations of the fluids 
and the interface between them. Beyond a certain threshold, the interface can form many 
kinds of standing wave patterns, including crystalline patterns and others which are more 



complex. This phenomenon was first studied by Faraday ( 1831 1 who noticed that the vi- 



bration frequency of the interface was half that of the forcing. The results of Faraday 
were confirmed by Rayleigh | ( 1883a|b I. Benjamin & Ursell (1954 1 carried out the first 
theoretical linear analysis of the Faraday waves, restricted to inviscid fluids. They decom- 
posed the fluid motion into normal modes of the container and showed that the evolution 
equation of each mode reduced to a Mathieu equation whose stability diagram is well 
known. 

In the 1990s, new behaviours of the interface were discovered, such as quasi-crystalline 
eight-fold patterns seen by Christiansen, Alstr0m & Levinsen (19921. By introducing 



a forcing which is the sum of two periodic functions with commensurable frequencies, 
Edwards & Fauve (19941 were able to produce twelve-fold quasi-patterns. Triangular 



patterns were observed by Miiller ( 1993 ) and superlattice patterns by Kudrolli, Pier 



& Gollub (1998), also using two- frequency forcing. Spatio-temporal chaos was studied 
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by Kudrolli & Gollub (19961, who also surveyed the occurrence of lattice patterns - 
stripes, squares or hexagons - as a function of viscosity and frequency. |Binks, Westra| 



& van de Water ( 1997 I demonstrated the dependence of the pattern on the depth of 



the layer. In addition to patterns or quasi-patterns, very localized circular waves called 
oscillons may occur, as seen by |Lioubashevski et al. ( 1996[ ). The Faraday instability is the 
first macroscopic system in which such structures have been observed. These discoveries 
endow the Faraday instability with a very great fundamental interest for understanding 
the natural formation of patterns. 

A number of theoretical or semi-numerical analyses were inspired by these experi- 
Kumar & Tuckerman ( 1994 1 extended the linear stability analysis oflBenjamin & 



mcnts. 



Ursell (1954) to viscous fluids. This analysis was experimentally confirmed by Bcchhoc- 



fer et al. ( 1995 ) and used by Kumar ( 1996 ) to predict cases in which the response would 



be harmonic rather than subharmonic. The method was extended by |Besson, Edwards fc| 



Tuckerman ( 1996[ ) to calculate the stability tongues in the case of two-frequency forcing. 



Integral equation formulations of the viscous linear stability problem were derived by 



Beyer & Friedrich (19951 and Miiller et al. (19971, who also studied the harmonic re- 



sponse case. Cerda & Tirapegui ( 1998 ) used the lubrication approximation and the WKB 
method to study shallow viscous layers, obtaining a Mathieu equation that was later used 
by Huepe et al. (2006) to derive analytic results about the response to multifrcquency 
forcing. 

Linear analysis provides no information about the shape of the patterns which ap- 
pear; other means are necessary to understand the occurence of a given pattern or the 
amplitude of stabilization. Weakly nonlinear approximations have been derived from the 
Navier-Stokes equations by Vinals and co-workers, e.g. Zhang & Vihals (19971, Chen 



& Vinals (19991 and by Skeldon & Guidoboni ( 2007[), foc using on the competition be 



tween different patterns. Vega, Knobloch & Martel ( 2001[ ) derived equations governing 
the interaction between Faraday waves and the mean flow. There has been a great deal 
of analysis of lattices, superlattices and quasi-patterns using equivariant dynamical sys- 
tems theory, as well as model equations designed to produce specific patterns, e.g. |Porter"^ 



Topaz & Silbcr ( 2004 ) . The approximation of quasipatterns in spatially periodic domains 



has also been addressed in Rucklidge & Silber ( 2009 ) . 

Investigation of the full nonlinear viscous problem, however, requires numerical sim- 
ulations, of which there have been very few up to now, specifically those of Chen fc| 



Wu 


(2000 


), 


Chen ( 


2002 


, |Murakami & Chikano 


(2001) 


Ubal, Giavedoni & Saita 


(2003 


) and 


O'Connor 


(- 


20081. 1 



(20081, all previous simulations have been two-dimensional. The most extensive sim- 
ulation thus far has been that of Chen & Wu ( 2000 1 and Chen ( 2002 1 , who used a 
finite-difference method applied to a boundary-fitted time-dependent coordinate system. 
At each timestep, the surface is advected and a new 2D grid, adapted to the surface, is 
recomputed. The amplitude of their numerically computed Faraday waves confirmed the 
weakly nonlinear analysis of Chen & Vihals (1999), including their prediction of a range 
of subcriticality. Their calculations also predicted qualitatively new phenomena, such as 
disconnected solution branches and slow modulated dynamics. 



Murakami & Chikano (2001 1 used a method similar to that of Chen & Wu (20001 and 



Chen (2002). Although they reproduced some features of the experiments by Liouba- 



shevski et al. (19961, their calculations were limited to accelerations only 0.5 % above 



critical. The investigation by Ubal et al. (2003) focused on the influence of liquid depth 
in two-dimensional simulations using a Galerkin finite-element method in transformed 
coordinates. In addition to comparing their linear stability predictions with those of |Ben-| 
iamin & Ursell (1954) andlKumar & Tuckermanl (1 19941), they calculated instantaneous 
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surface profiles and velocity fields, as well as the temporal evolution and spectrum. Valha 
et al. (2002) examined the response of a liquid layer in a vertical cylindrical vessel us- 
ing the MAC (Marker-and-Cell) method of Harlow & Welch (19651. Surface tension was 
treated by the continuum surface force model of iBrackbill, Kothe & Zemach (19921. 



O'Connor] ( 2008 ) conducted numerical simulations using an ALE (Arbitrary Lagrangian- 



Eulerian) spectral-element code in both two and three dimensions; a visualization of a 
square pattern was presented. 

The hexagonal patterns, quasi-patterns and oscillons which motivate our investiga- 
tion are intrinsically three-dimensional, and have never been calculated numerically from 
the fluid-dynamical equations. Here we report on the results of fully nonlinear, three- 
dimensional simulations of Faraday waves using a finite-difference front-tracking method. 
In the classic Faraday problem the lighter fluid is usually taken to be air whose effects 
can be neglected. However, in contrast to the previously cited investigations, the numer- 
ical method described here solves the Navier-Stokes equations for the general case of 
two distinct superposed fluids. The capability of the method to simulate the motion of 
both fluids is important in that it permits comparison of numerical results with those 
of certain experimental configurations, namely those of Kityk et al. ( |2005[ ) where the 



lighter fluid cannot be ignored. These experiments were the first to provide quantitative 
measurements of the complete spatio-temporal Fourier spectrum of Faraday waves and 
thus form an excellent basis for quantitative comparison with our numerical results in 
the nonlinear finite amplitude regime, i.e. interfaces with steep slopes. 

In the next two sections of this article, we present the hydrodynamic equations that 
govern the Faraday instability and then describe the computational method. The two 
sections following these are dedicated to the comparison of our results with the linear 
theory of Kumar & Tuckerman (1994) and with the experiments of Kityk et al. (2005, 
2009). After comparing numerical and experimental spatio-temporal spectra for squares 
and hexagons, we present the three-dimensional velocity field for the hexagonal pattern. 



2. Equations of motion 

The mathematical model of the Faraday experiment consists of two incompressible 
and immiscible viscous fluids in a three-dimensional domain x = (x,y,z) G 5R 2 x [0, h], 
bounded at z = and z = h by flat walls. The two fluids, each uniform and of densities 
pi, p2 and viscosities p\ and p2, initially form two superposed horizontal layers with an 
interface between them. This two-dimensional interface is defined by x' = (a;, y, £(x, y, t)). 
Within the parameter range we wish to simulate, the height £ remains a single-valued 
function of (x, y, t). 

The container is shaken vertically in z. In the reference frame of the container the 
boundary conditions for the fluid velocities u = (u,v,w) are 

u(x,y,0,t) = 0, (2.1a) 
u(x,y,h,t) = 0. (2.16) 

The gravitational acceleration g is augmented by a temporally periodic inertial acceler- 
ation 

G = (acos(wi) - g)e z , (2.2) 

where a is the amplitude of the forcing and to is its frequency. 

The Navier-Stokes equations for incompressible, Newtonian fluids are 

= -Vp+pG+V ■ /i(Vu+ Vu T ) +s, (2.3a) 
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V-u = 0. (2.36) 
Here p is the pressure and s is the capillary force (per unit volume) and is defined below. 



Equations (2.3a) and (2.36) are valid for the entire domain, including the interface, in 
spite of the fact that the density and viscosity change discontinuously and the surface 
tension acts only at the interface. In this single-fluid formulation, the density and viscosity 
fields are defined in terms of the densities and viscosities of the two fluids 

P = Pi + {P2 - Px)H, (2.4a) 
H = fit + (M2 - Hi)H, (2.46) 
with the aid of a Hcavisidc function, 

H{x-*) = \°**Z% X ' V ' t 2 , (2.5) 
]_ 1 if z ^ C{x,y,t) 

where we recall that x = (x, y, z) is a point anywhere in the three-dimensional volume 
and x' = (x,y, C,{x,y, t)) is the vertical projection of x onto the interface. The capillary 
force is 

s = / (TKnS(x-x')dS, (2.6) 

JS'(t) 

where a is the surface tension coefficient, assumed to be constant, n is the unit normal 
to the interface (directed into the upper fluid) and k its curvature. 5 (x — x') is a three- 
dimensional Dirac distribution that is nonzero only where x = x'. S'(t) is the surface 
defined by the instantaneous position of the interface. 

To complete the system of equations we need an expression for the motion of the 
interface. One such expression can be easily derived by noting that mass conservation in 
an incompressible flow requires Dp/Dt = 0, which in view of the discontinuous density 



field (2.4a I, is equivalent to 

DH/Dt = 0. (2.7) 

Thus the interface is represented implicitly by H and advected by material motion of the 
fluid. 



3. Computational methods 

The computational domain is a rectangular parallelepiped, horizontally periodic in x 
and y and bounded in z by flat walls for which we impose no-slip boundary conditions. 
The entire domain is discretized by a uniform fixed three-dimensional finite-difference 



mesh. This mesh has a standard staggered MAC cell arrangement (Harlow & Welch 



1965 1 where the u, v and w velocity nodes are located on the corresponding cell faces 



and scalar variables are located at the cell centres. Each cell is of dimension Aa; x Ay x Az. 

Within the domain, the two distinct immiscible fluids are separated by a two-dimensional 
interface which is discretized by a second mesh as sketched in figure [T] This moving and 
deformable mesh is composed of triangular elements whose motion is treated by a front- 
tracking/immersed-boundary method ( |Peskin|1977 Tryggvason et aL|2001[ ). Because we 
have assumed that ((x, y, t) is single valued, the nodes of the mesh can be fixed in x and 
y and only their vertical displacements need to be calculated, which is a considerable 
simplification to the general front-tracking method. 

After setting appropriate initial and boundary conditions, the computational solution 
algorithm for each timestep is composed of three main phases. First, the interface is 
advected and the density and viscosity fields updated according to the new interface 
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Figure 1. Spatial discretization of the domain. 



position. The capillary force s is then calculated. Finally, the velocity and pressure are 
found by means of a standard projection method. Each of these steps is described below. 

3.1. Advection of the interface 
Purely Eulerian interface methods such as Volume of Fluid (Hirt & Nichols 1981) or 



Level Set (Osher & Sethian 1988) use a form of (2.7 1 to advect a scalar field such as 



H, or a level-set function that implicitly represents the interface. However, in the front- 
tracking approach that we use here, the interface markers themselves (the nodes of the 
triangular mesh) are advected. H is then constructed from the position and geometry of 



the interface. Taking (2.7 1 as a starting point, we develop an equivalent expression for 



the vertical displacement of the triangular interface mesh. 



The material derivative of (2.5) gives 



BH 



VH 



Dx 



V'H ■ 



Dx' 

DP 



where V' = <9x< and 



V/f = -V'H = [ n<J(x-x')dS. 

JS'{t) 



(For a derivation of (3.2 1 see Tryggvason et al. 2001.) Factoring (3.1l by ViJ 



VH 
~Dt 



VH- u 



Dt 



(3-1) 



(3.2) 



(3.3) 



The right-hand side of (3.3 1 can only be zero everywhere, including on the surface, x = x', 
if 

Dx' , , s , N 

^ = u(x',i), (3.4) 
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which establishes the material motion of the explicit interface representation x'. Further- 
more 

| = (I-Vx').u(x',i), (3.5) 

where I is the identity tensor. The specific choice of x' = (x, y, £(x, y, t)) made here gives 
for Vx' 

1 
1 

d( d{ 



dx dy 







(3.6) 



With this, (3.5 1 leads to the specific displacement relations: 

dx/dt = 0, 



(3.7a) 



dy/dt = 0, 
d( 



w 



dC 



(3.7b) 
(3.7c) 



dC 

dt dx " dy 

The application of this advection to the triangular mesh we use for tracking the interface 
is straightforward. At each vertex the horizontal displacement is zero and for the vertical 



displacement we compute a first-order approximation to (3.7c 



c 



ra+1 



At 



dx 



dy [ e) - 



(3.8) 



The superscripts n and n + 1 denote, respectively, the old and new time levels. The 
derivatives on the right-hand side are evaluated using a simple upwind scheme which 
requires the usual CFL (Courant-Friedrichs-Lewy) time step restriction. The vertical 
displacement of the interface mesh requires knowledge of the velocities at the element 
nodes x' a . These in general do not coincide with the Eulerian grid nodes x^-fe, whose 
indices correspond to discretized coordinates along the respective directions x, y and z. 
The problem of communicating Eulerian grid velocities to the clement nodes is over- 
come by interpolation between the two grids as is typically done in front-tracking and 
immersed-boundary methods. Here we use the particular interpolation 



u(x' e ) 



u(x ijk )5 h (x ijk - x'JAxAyAz. 



(3.9) 



The kernel 5^ is a smoothed version of the three-dimensional Dirac delta function with 
compact support of four grid nodes in each direction (for details of the front-tracking 
method, see Tryggvason et al. 2001, and for the immersed-boundary method, see Peskin 



1977). In (3.9) the weighted information collected from nearby Eulerian grid nodes is 
interpolated to a given element node. 



We now seek to update the density and viscosity fields needed in (2.4 1, which require 



H. The equation for H , based on the updated values of x' and n, is formulated by taking 



the divergence of (3.2 1: 



l H = V- [ nS(x-x')dS, 
Js'(t) 



IS'(t) 
H(x,y,0,i)=0, 



(3.10a) 
(3.106) 



H(x,y,h,t) = 1. 



(3.10c) 
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Figure 2. A trian gular element of the interface mesh illustrating the action of the capillary 
forces according to (3.131. For the shaded triangle the forces act perpendicular to the triangle's 
edges (unlabelled solid arrows), the dashed arrows are the corresponding forces on the edges 
of neighboring triangles. The net capillary force at the shared edge between any two triangles 
(the sum of the solid and dashed vectors) is directed into the fluid on the concave side of the 
interface. 



The discretized version of the Poisson problem (3. 10a I is 



(3.11) 



where standard central differencing is used for the gradient and divergence operators. 
This numerically calculated Heaviside function is a smoothed transition from to 1 
across a distance of 4 grid cells in the direction normal to the interface. In contrast to 
(3.9 1, the summation above serves to distribute weighted information from an clement 



node to nearby Eulerian grid nodes. Since an clement is triangular, its vertices lie in 
the same plane, its normal vector is unique and the three tangent vectors are simple to 
calculate: 

*2Xti 

n c = - - (3.12) 

||t 2 x till v ; 

where tx and t 2 are the tangents on two distinct edges of the triangle (see the sketch in 
figure [2|. W e solve (3.11) by fast Fourier transform. Finally, p n+1 and /i n+1 are updated 
using \2Aa\ and \2Ab\. 



3.2. Capillary force 

From (2.6), the capillary force involves the curvature of the interface and its normal 



vector. However, from a computational point of view, curvature is a difficult quantity 
to compute accurately. It is more accurate and physically appealing to calculate the 
force pulling on the edge of each individual triangular surface element and then sum 
the contributions for all the elements over the surface. For a given surface element e of 
surface area S A and perimeter SI, we can write: 



s e = (7 / Kn dA, 

J 8 A 



a (n x V) x n dA, 

JSA 



(3.13) 
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where the last integral represents the sum of the capillary forces exerted around the 
element perimeter. As sketched in figure [2] the directions of these forces are oriented 
along the surface and normal to the element's edges. The net force at the shared edge 
between any two triangles (the sum of the solid and dashed vectors) is directed into the 
fluid on the concave side of the interface. Following Peskin's immersed boundary method 
( Peskin|[l977 ), the discrete version of ( |2.6[ ) becomes 



ft = ^ S e^h(^i]k - x' e ), 



(3.14) 



where we use the same smoothed 6h as in (3.9 1 and (3.111. Thereby, several interfacial 



elements contribute to the calculation of the force applied to a single Eulerian node, and 
a single element influences more than one Eulerian node. 

3.3. Solution of the Navier-Stokes equations 

The Navier-Stokes equations are solved by a projection method (Chorin 1968; Temam 
1968) with incremental pressure correction (Goda 1979) applied to a finite-difference 
scheme which is first order in time and second order in space. In addition a semi- implicit 
scheme is chosen for the velocities to relax the stability restriction on the time step due 
to viscous diffusion. All spatial derivative operators are evaluated using standard centred 
differences, except in the nonlinear term where we use a second-order ENO (Essentially- 
Non-Oscillatory) scheme ( |Shu fc Osher 1989 Sussman et al. 1998[ ). (For an overview 
of projection methods for the incompressible Navier-Stokes equations, see Guermond, 
Minev & Shen 2006.) The time stepping algorithm is thus 



,n+l 



At 



u n • Vu" + — - V • u" +1 (Vu + Vu T 



n+1 



1 



with the boundary conditions on the top and bottom walls 

u" +1 L = 0. 



r.n+1 



G 



n+1 



(3.15) 



(3.16) 



In (3.151, p, p and s depend on x via (2.4-2.6) and have already been updated by (3.8 1. 



We decompose the solution of (3.15) in three steps. The first step is a semi-implicit 



calculation of an intermediate unprojected velocity u, involving only velocities and their 
gradients: 

1 



u 



At 



= — u 



Vu"+^ I V. /J " +1 (Vu+Vu T ) 



P 



(3.17a) 



u| r = 



(3.176) 

In the second step, we include the capillary, acceleration and old pressure gradient terms 
to calculate the unprojected velocity u*: 

U 11 Qn+l b 1 



At ~ ' p n+1 p" 
Finally we perform a projection step to find the divergence free velocity u™ +1 : 



1+1 ^ 



(3.18) 
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n+l 



= o, 



u" +1 -n| r = 0. 



Equations (3.19) imply the following elliptic problem for the pressure increment 



At 



V • 



1 



n+l 



P 



dip 



n+l 



dn 



= 0, 
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(3.196) 
(3.19c) 

(3.20a) 
(3.206) 



which we solve with an iterative biconjugate gradient stabilized algorithm (Saad 1996). 
In the horizontal directions, periodic boundary conditions are imposed on the velocity 
and pressure, thus excluding a net horizontal pressure gradient. For the simulations we 
will present here, this choice is consistent with the requirement of no mean horizontal 
flux in a large but bounded container. 



We note that in the implicit solution of (3.171, we apply the same biconjugate gradient 
stabilized solver used for the pressure to each component of u = (u, v, w) separately. 
Thus only the diagonal terms of the diffusion operator are treated fully implicitly. The 
off-diagonal terms are treated quasi-implicitly in that the newest available values of 
(u,v,w) are used in the evaluation of the cross derivatives. To ensure symmetry, we 
permute the order of solution for each component. 



4. Results: linear analysis 

4.1. Floquet analysis 

In the absence of lateral boundaries, the equations are homogeneous in the horizontal 
coordinates and the solutions can be represented by a spatial Fourier transform: 

C(x, y,t) = C(M)c ik ' x (4.1) 

k 



The linear instability of the interface between two fluids is described by ( 2.2 1 — (2.7 1 lin- 
earized about a zero velocity field and flat interface £ = (£). The linearized equations 
depend only on the wavenumber k = ||k|| of each wave and not on its orientation and 
hence the coefficient of e lkx can be written as ((k,t); additionally the dynamics of each 
t) is decoupled from the others. Linear partial differential equations with constant co- 
efficients have solutions which are exponential or trigonometric in time. For the Faraday 
instability t) is instead governed by a system of linear partial differential equations 
with time-periodic coefficients, i.e. a Floquet problem, whose solutions are of the form 

C(fc, t) = e^ +iQ ^7(fc, * mod T), (4.2) 

where T = 2n/u>, 7 is real and a £ [0, 1[. The Floquet modes, 

00 

f(k,t mod T)= J2 /n(*)e inwl , (4-3) 

n=— 00 



are not trigonometric, but remain periodic with fundamental frequency lj. Thus, the 
linearized behavior for a single mode is 

00 

C(x,y,t) = e ik - x e (7+imj) * J] f„(k)e inbJt . (4.4) 
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Analogous expressions hold for the velocity u. 



Equation (4.2) shows that if 7 is non-zero or a is irrational, the evolution of the 
interface motion is not periodic. A non-zero 7 indicates that the motion grows or decays 
according to the sign of 7. An irrational a yields a quasi-periodic evolution function. 
For the Faraday instability, it can be shown that a can only take two values: and 1/2. 
As the imposed acceleration a is increased, one encounters regions in the (fc, a) plane in 
which 7 > for one dominant temporal frequency, jw/2, where j = 1,2,3, ... (see figure 
I3J. Within each instability tongue, the amplitude of the mode grows exponentially. These 
tongues are called harmonic if a = Q and subharmonic if a — 1/2. As k is increased, 
one encounters an alternating sequence of subharmonic and harmonic tongues, which are 
bounded by neutral curves (k, a c {k)) on which 7 = 0. On the neutral curves, the solutions 
are periodic: 

00 

C(fc,i)= fn(k)e^ n+ ^ ut , subharmonic case, (4.5a) 

n— — OO 

OO 

C(M)= fn(k)e [nuJt , harmonic case. (4.56) 



4.2. Computation of the neutral curves 
We first compare our numerically calculated instability thresholds with those found by 



Kumar & Tuckerman (1994) for the same parameter values. The physical parameters 
are p\ = 519.933 kg m -3 and \i\ = 3.908 x 10~ 5 Pas for the lower fluid and P2 = 
415.667 kgm~ 3 and [12 — 3.124 x 10 -5 Pas for the upper fluid. The other parameters 
are a = 2.181 x 10 _6 Nm and g = 9.8066 ms~ 2 . The frequency of the forcing is 
u/2n = 100 Hz and thus its period is T = 0.01s. The capillary length is defined as 
l c = y 'o /dp! — P2\g)- The container height is taken to be 5l c = 0.231 mm, and the 
interface, when unperturbed, is equidistant from the top and bottom boundaries. We 
consider several wavenumbers k and set the x dimension of the box in each case to one 
expected wavelength A = 2n/k, i.e. to between 0.074 and 0.224 mm, as listed in Table 
[l] We can estimate the importance of various physical effects for these parameters by 
defining dimcnsionless quantities with length and forcing period T. The Bond number 
Bo = (kl c )~ 2 = \pi — P2I g/i^k 2 ) measures the relative importance of gravitational to 
capillary effects and ranges between 0.0649 and 0.598. The Reynolds number Re — 
p/(pk 2 T) is a nondimensional measure of viscous damping and ranges between 0.184 
and 1.70 for both fluids. 

We have computed the critical acceleration from our fluid-dynamical simulation for 
the wavenumbers listed in Table [l] Initially, the interface is sinusoidal with wavevector 
k parallel to the £-axis and the velocity is zero. Moreover, to ensure that the solution 
corresponds initially to the linear solution, we require the amplitude of the interface 
displacement to be small compared to A. In order to maintain a roughly cubic mesh and 
a minimum x-resolution of about 50 grid cells per wavelength, we vary the resolution 
in the z direction between 126 and 144 points. Since k points along the x direction, £ 
does not depend on y (neither do the velocity nor the pressure) and so the size of the 
domain and resolution in y are arbitrary. The acceleration a is taken near a c {k), the 
expected critical acceleration corresponding to each wavenumber. At the threshold, the 
flow undergoes a pitchfork bifurcation. Since the growth rate is proportional to a — a c 
close to the neutral curve, it is sufficient to find the growth rates for two values of a and 
to interpolate linearly between these points. 

In figure [3] we plot the values of a c obtained from our fluid-dynamical simulation 
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k (mm J ) 



Figure 3. Critical acceleration a c /g as a function of the wavenum ber k. The solid curves 
represent the neutral curves obtained by |Kumar fc Tuckermari] ( |1994| . The a c found with the 
simulation are indicated by the circles. 



k( mm ) A( mm) No. of gridpoints in x/z a c /g (Theor.) a c /g (Comp.) Error(%) 

28 0.224 124/128 

32.5 0.193 96/128 

35 0.180 100/128 



0.131 72/126 
60.9 0.103 56/126 
85 0.074 48/144 

Table 1. Comparison of the computed a c with Floquet theory for various wavenumbers k. 
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for several values of k, along with the curves (k,a c (k)) obtained from the method of 
Kumar & Tuckerman (19941. Figure [3]shows that these thresholds are in good agreement, 



despite whatever inaccuracies in a c are introduced by spatial discretization and linear 
interpolation. The relative error in the critical acceleration at the conditions previously 
stated is of the order of a few per cent as shown in Table [T] 

The results suggest that there is a A: below which calculation of the growth rates is 
not possible. Some zones of the diagram are not accessible because a domain of width 
2n/k necessarily accommodates all wavenumbers which are integer multiples of k up 
to the resolution limit n/kAx. The coefficients of the Fourier expansion of the initial 
condition differ slightly from zero due to finite-difference spatial approximations and if 
the growth rate of one of these is greater than that of k itself, then it will quickly come 
to dominate k. This difficulty is exacerbated by the fact that several forcing periods 
are required for 7 to stabilize. Then the amplitude, whose evolution was expected to 
be almost periodic, starts to rise before the precise determination of 7 is possible, for 
example in the range of k between and roughly 15 mm" 1 . As we see in figure U the 
critical forcing is substantially lower for one of its multiples closer to 32.5 mm" 1 . The 



12 



N. Perinet, D. Juric and L. S. Tuckerman 



amplitude corresponding to this wavelength, although initially negligible, increases and 
rapidly dominates the mode we wish to study, making the calculation of 7 unfeasible. 
In contrast, for k — 48 mm , the growth rate did not vary significantly after having 
reached a value near zero (relative fluctuations of about 0.1% of the growth rate's limit 
value were recorded after the stabilization). 

4.3. Temporal profile of a mode 



We recall from section 4.1 that the time dependence of a Floquet mode is not sinusoidal. 
As a further validation, we can compare the results of our fluid-dynamical simulations 
to the entire temporal behavior over a period. This is a stronger validation than merely 
predicting the threshold since it provides a comparison at every time instead of once per 
period. 

In figures |ij|6| we plot the deviation £ — (Q from the flat interface as a function of 
time at a fixed spatial location from our fluid-dynamical simulation, for values k = 48, 
60.9 and 85 mm^ 1 belonging to the first three tongues. On the same figures, we plot 



the behavior of (4.5 1, where the temporal coefficients f n (k) of the Floquet modes have 
been calculated by the method in Kumar & Tuckerman ( |1994 1 . The value of a is set to 



the interpolated critical acceleration a c {k), so the oscillations approximately retain their 
initial amplitude as long as they remain small. The comparisons in figures 4 r (> show a 
nearly perfect agreement. The differences observed initially, due to the phase difference 
between the initial conditions and acceleration, vanish remarkably quickly, in well under 
one period of oscillation of the container. 

Figures [4jj6] correspond to tongues jui/2, with j =1, 2, 3, respectively, which show j 
zero crossings per forcing period T. Odd (even) values of j correspond to subharmonic 
(harmonic) oscillations, with period 2T (T). The temporal spectrum f n (k) becomes richer 
as k increases, leading to increasingly more complex modes, as can be observed by com- 
paring figures [4||6] This strong anharmonicity of the curves is due to the increasing 
contribution of higher frequency trigonometric functions to the Floquet modes as a in- 
creases. The Floquet mode corresponding to k c = 32.5 mm -1 , with the smallest value of 
a = a c , should be closer to trigonometric, with a fundamental frequency of u/2. 



5. Results: nonlinear analysis 

In the full nonlinear evolution of the interface for a > a r 



the amplitude of the interface 



height grows in time until nonlinear terms in (2.2 1— (2.7 ) become important. After that, 
the mode whose linear growth rate is maximal gives rise, via nonlinear resonances, to a 
series of other discrete modes, selected according to the magnitudes and orientations of 
their wavevector k. This selection is responsible for the formation of patterns that will 
be the object of our further validations. We seek to compare our calculations with the 
experimental results of Kityk et al. (2005, 2009) where quantitative data concerning the 
Fourier spectrum C(k, f) are available for squares and hexagons. 



We run our numerical simulations with the same experimental parameters as Ki- 

12 Hz (T = 0.0833 s), p x = 1346 kgm~ 3 , p x = 7.2 mPas 
949 kgm -3 , p 2 = 20 mPas for the upper fluid. The surface 
= 35 mNm -1 , the total height of the vessel is 1.0 cm and the 



tyk et al. (20051: co/2n 



for the lower fluid and p 2 - 
tension at the interface is a 
mean height of the interface, the initial fill height of the heavy fluid, is (Q = 1.6 mm (with 
some uncertainty; see below) . The Floquet analysis for these parameters yields a critical 
wavelength of A c = 2n/k c = 13.2 m m and a critical acceleration of a c = 25.8 ms -2 . Here, 
the Bond number defined in section 4.2 is Bo — \pi — p 2 \ g/(ak^) = 0.49. The Reynolds 




Figure 4. Linear evolution of the surface height deviation £(t) — (£) for k — 48 mm -1 , in 
the first instability tongue. Our simulation results are plotted with symbols and those derived 
from a Floquet analysis with the solid line. The height and time are nondimensionalized by the 
wavelength A = 2n/k and forcing period T, respectively. 

0.006 
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0.004 



0.002 



-0.002 




Figure 5. Linear evolution of the surface height deviation — (£) for k — 60.9 mm 1 , in the 
second instability tongue. Same conventions used as figure [4] 
0.03 r 




-0.03 



Figure 6. Linear evolution of the surface height deviation ((t) — {() for k = 85 mm , in the 
third instability tongue. Same conventions used as figure [4] 
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number Re = p/(fj,k^T) is Re\ — 9.9 and Re^ — 2.52 for the lower and upper fluid, 
respectively. 

Rather than starting from a sinusoidal interface, we chose to add two-dimensional 
white noise of small amplitude to (0 to define the initial interface height £(x, y, t = 0) 
in order to excite every mode allowed by the box's horizontal dimensions and number of 
cells. It is thus possible to check that the correct critical mode (that whose growth rate 
is maximal) emerges from the linear dynamics. In order to reproduce the experimental 
results in a computational domain of a minimal size, the dimensions in x and y of the 
box must correspond to the periodicity and symmetries of the expected pattern. The 
minimal required resolution along these directions has been found to be between 40 and 
50 cells per wavelength. The number of triangles used to represent the interface is 16 
times the total number of horizontal gridpoints. The number of cells in the z direction 
is taken so that imn((x, y,t) is greater than about the height of 3-5 cells. The required 

vertical resolution thus varies with the forcing amplitude. The initial velocity is taken to 
be zero. 

5.1. Square patterns 



To compare with the experiment of Kityk et at. ( 2005 1 for their square patterns, we 
choose the same forcing acceleration, a = 30.0 ms - ^. Our box has horizontal dimensions 
which we take both equal to 2n/k c . The timestep is At = 2.78 x 10~ 4 s. Figures [7] and [5] 
represent examples of the patterns obtained at saturation under these conditions and are 
taken from the same simulation at the two instants shown by the two arrows in figure 



10 The symmetries characterizing the squares (reflections and n/2 rotation invariance) 



are clear, showing a first qualitative agreement with Kityk et al. (20051 where both 
structures were observed. The pattern oscillates subharmonically, at 2T, where T is the 
forcing period. Figure [7] is taken when the interface attains its maximum height, while 
figure [8] is taken at a time 0.24 x 2T later. At this later time, we observe the dominance 
of a higher wavenumber, which will be discussed below. 

Further quantitative investigations of the patterns involve the spatial Fourier transform 
of the interface height. In the case of square patterns, the distribution of the spatial modes 
is shown in figure[9] The modes with non-negligible amplitude are ztk c e x and ±k c e y , with 
|k| = k c and amplitude A(k c ); ±2k c e x and ±2k c e yi with |k| = 2k c and amplitude A(2k c ); 
and fc c (±e x ± e y ), with |k| = y/2k c and amplitude A(y/2k c ). (For a square pattern, the 
amplitude of each mode is identical to that of each of its images through rotation by any 
integer multiple of tt/2.) The interface height is written as: 

4 4 

C(x,f)= (C> + A(k c ,t) 5>^- x + A(2fe c ,t) ]Te i2fe ^- x 

4 

+ A(V2k c , t) e iv/5 M x + higher order terms, (5.1) 

where ej = e x cos(irj/2) +e y sin(7rj'/2) and e^- = e x cos(7r/4 + irj/2) +e y sin(7r/4+ nj/2) 
for j = 1, . . .4. We have chosen this notation, rather than C(^i^) as used in equation 



(4.1|, to facilitate comparison with Kityk et al. (2005, 2009). 

We have compared the evolution of the three principal spatial modes (figure 10 1 



and their temporal Fourier transform (figure 1 1 1 with the experimental results (Kityk 



et al. 2005). Here we turn the reader's attention to the recent erratum by Kityk et al. 



(2009) for correct quantitative comparisons of the spectra. 



Figure [10| compares the experimental evolution of each spatial wavenumber to numcr- 
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Figure 7. Example of square pattern. Height of interface as a function of the horizontal coordi- 
nates, at the instant corresponding to first arrow of figure [To| when height is maximal. Resolution 
in x, y, z directions: 80 x 80 x 160. Note that the vertical scale is stretched with respect to the 
horizontal scale. Each horizontal direction in the figure is twice that of the calculation domain. 



4-v 
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Figure 8. Example of square pattern, at the instant corresponding to second arrow of figure 

[TO] time 0.24 x 2T after figure [7) 




ical calculations for two different mean heights (£) = 1.6 mm and (£) = 1.7 mm. Our 
calculations show that the results depend strongly on ((), which is the initial fill height 
of the heavy fluid. Our discussions with Kityk and Wagner (A. Kityk & C. Wagner, 
private communication) indicate that this is true as well in the experiments, and also 
that 0.1 mm is within the experimental uncertainty for their mean height. Thus we chose 
to vary (£), in preference to other parameters, in order to check whether the range of 
amplitudes caused by experimental uncertainties includes those obtained numerically 

The main features found in Kityk et al. ( 2009[ ) are recognized in figurejlOj In particular, 
both the fundamental periodicity of each mode (harmonic or subharmonic) and the form 
of each numerical curve in figure [10] are very similar to the experimental data. The 
amplitudes and the phases are also quite close. Most of the experimental amplitudes are 
bracketed by the numerical ones. Thus, they lie in the interval of amplitudes allowed 
by the range of uncertainties which is surely underestimated since only the uncertainty 
in (£) has been taken into account. A(k c ) crosses zero at times different from the two 
higher wavenumbers, A(2k c ) and A(\^2k c ). At these instants, the higher wavenumbers 
dominate the pattern. In particular, the pattern of figure [8] taken near the second arrow 
in figure 10 when A(k c ) is low, contains more peaks than that of figure [7] taken when 
A(k c ) is high. The large ratio between the amplitude of k c and the others makes this 
phenomenon very short-lived. 

Figure [TT] shows the temporal Fourier decomposition of the curves in figure |10| These 
spectra for the experiment Kityk et al. ( 2009 1 and for the computation are quite similar 
too. All of the square patterns that we have observed, once saturation is attained, remain 
so for the entire duration of the calculation. 



We present a brief numerical grid convergence study in figure 12 All qualitative fea- 
tures, such as the square symmetry, were observed with each of the three resolutions 
chosen, despite the coarseness of the 20 x 20 x 40 and 40 x 40 x 80 grids. With increasing 
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10 10.05 

time(s) 

Figure 10. Temporal evolution of the amplitudes of the spatial modes with wayenum bers k c , 



2k c and yf2k c . Solid curves represent the experimental results of Kityk et al. (20091, dashed 
curves and crosses represent numerical results for (£) = 1.6 mm and (Q = 1.7 mm, respectively. 
Resolution in x, y, z directions: 80 x 80 x 160. Arrows, from left to right, show the time at which 
figures [7] and [8] have been plotted. 
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Figure 11. Tem poral Fourier trans form of the amplitudes in figure [To| Circles indicate experi- 
mental results of |Kityk et al. (20091, while crosses and plus signs indicate numerical data with 
(0 = 1.6 mm and (Q = 1.7 mm respectively. 
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-0.2 b i i i i i J 

27.4 27.45 27.5 27.55 27.6 

time(s) 

Figure 12. Temporal evolution of the amplitudes of the spatial modes with wavenumbers k c , 2k c 
and V2k c for square patterns. Study of the convergence with three different spatial resolutions. 
Circles indicate a resolution (in x, y, z directions) of 20 x 20 x 40, dashed curves 40 x 40 x 80 and 
continuous curves 80 x 80 x 160. The timestep is the same, At = 2.78 x 10 -4 s, for all curves 
shown. 



resolution, the principal spatial modes converge to the experimental curves shown in fig- 
ure [lO] with only a small difference between the curves with the two highest resolutions. 
The order of numerical convergence of the maximum and minimum of the amplitudes 
of each of the three modes in figure [12] shows that the convergence is between first and 
second order, which is expected to be the case with the immersed-boundary method. In 
particular, we would expect that a further doubling of the resolution would change the 
results by at most 4 % for the principal k c mode. 



5.2. Hexagonal patterns 
When the amplitude of the forcing acceleration a is further increased, the modes can 



reorganize. The symmetries change and, in the experiments of |Kityk et al. (20051, the 
initial square pattern becomes hexagonal. Though k c remains constant, the horizontal 
dimensions of the minimal computational box necessary to support the periodic pattern 
must change too. These dimensions become An/k c in y and 47t/(v / 3/c c ) in x, as shown 



in figure 13 The wavevector lattice for hexagonal patterns is shown in figure 14 The 
principal modes are again of three amplitudes: k c , 2k c and %/3fc c . When a pattern is 
hexagonal, a mode will have the same amplitude and temporal behaviour as each of its 
images through rotations by any integer multiple of 7t/3. The interface height is thus 



C(x, t) = (C) + A(k c , t) " x + M^kc 

3=1 



6 
3=1 



2k c 
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Figure 13. Boxes supporting the periodic patterns in the square and hexagonal cases. In 
black, the borders of the box. Light lines, pattern contained by each box; A = 2n/k c . 



A(V3k c ,t) 



3 = 1 



higher order terms, 



(5.2) 



-e y sin(7rj/3) and e^- = cos(7r/6 + 7rj73) +e y sm(7r/6 + 7rj/3). 



where e., = e x cos(nj/3) - 
for j = 1, .. .,6. 

Our simulations are carried out at acceleration a = 38.0 ms -2 and mean height (Q = 
1.6 mm. We have used two different initial conditions: a rectangular pattern, and also 
white noise, as in our previous simulations of the square patterns. In both cases, hexagons 
emerge and saturate. The results shown below are those that emerge from the white noise. 
The time step varies during the calculation, depending on the viscous diffusion limit and 
the CFL. The spatial resolution is 58x100x180 in the x, y, z, directions, respectively. 

In figures |T5f[T8| we show visualizations of the patterns at four instances in time. A 
movie of the temporal evolution of the hexagon pattern over one subharmonic oscillation 
is available in the online version of this article. The 7t/3 rotational symmetry confirms that 
the rectangular numerical grid does not forbid the formation of hexagonal patterns, which 
are not aligned with this grid. The patterns reproduce several prominent features from 
the visual observations of hexagons in the experiments. For example, one can observe 



in figures 15 and 17 the up and down hexagons shown in the experimental snapshots 



(figure 10 of Kityk et al. 2005). The pattern in figure 18 when the surface elevation is 



minimal, is dominated by wavenumbers higher than k c , as is also the case in figure 10 



of |Kityk et al. (20051. This is reflected by the disappearance of A{k c ) and the resulting 
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Figure 14. Lattice formed by the spatial modes comprising a hexagonal pattern. The principal 
modes, with wavenumbers k c , 2k c and y/3k c , involved in later quantitative investigations are 
indicated by hollow black circles. The labelled triangles illustrate resonance mechanisms leading 
to harmonic contributions to higher wavenumbers. 



dominance of A(2k c ) and A(^/3k c ) at the corresponding instant in the spectral timeseries 
of figure |19| This apparent wavenumber increase is analogous to that which occurs for 
the squares, shown in figures |8] and [10] 



The spectra from experiments and simulations are represented in figures 19 and 20 



Given that experimental uncertainties concerning the hexagons are greater than for the 
squares (A. Kityk & C. Wagner, private communication), the agreement is remarkable. 
The principal mode is well reproduced while the other two modes show rough agreement . 
It is striking that, in contrast to square patterns, every wavevector is a superposition of 
harmonic and subharmonic temporal modes, so that each has temporal period 2T. This 
phenomenon was explained by Kityk et al. ( 2005 1 as a spatio-temporal resonance as 
follows. In the case of the square lattice, two critical subharmonic modes (e.g. k c e x and 
k c e y ) interact to yield a higher wavenumber harmonic mode (e.g. k c (e x + e^)). In the 
hexagonal case, two critical subharmonic modes (e.g. —k c e y and fc c (V3e x — e y ) /2) interact 
to y ield a higher wavenumber harmonic mode (k c (^/3e x — 3e y )/2), as in triangle I of figure 
Il4l Further interaction of this mode with a critical subharmonic mode (k c (V3e x + e y )/2) 
yields subharmonic contributions to the higher spatial wavenumber mode (k c (y/3e x — e y )), 
as shown in triangle II. Other quadratic interactions between critical subharmonic modes 
can contribute to a third harmonic mode of wavenumber k c (triangle III) . 

In addition to the interface height, our simulations also produce the entire velocity field, 
which is the focus of figures [2lf|23| These figures show the velocity fields on horizontal 
planes at three instants spanning the oscillation period of a hexagonal pattern, as well 
as the vertical velocity on the interface. Figures 21 22 and 23 correspond approximately 
to the visualizations of figures 15 [16] and [18] where the structures are more visible since 
the interface has been repeated periodically in the horizontal directions for clarity. The 
parameters are the same as those given previously, except that the acceleration a has 
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Figure 15. Snapshot of hexagonal pattern, taken when height of the interface peaks is 
maximal. Each horizontal direction is twice that of the calculation domain. 
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Figure 16. Snapshot of hexagonal pattern taken t — 0.3 x 2T after instant of maximal 

interface height. 



been decreased to 36.0 ms~ 2 , and the number of triangles used to represent the interface 
has been increased to 64 times the total number of horizontal gridpoints. 

Figure 21 is taken at t = 0.07 x (2T), just after the interface reaches its maximum height 
(at t = 0), when the peaks are beginning to descend. Consequently, the fluid converges 
horizontally towards the interface peaks, then descends dramatically below them. The 
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Figure 17. Snapshot of hexagonal pattern taken t = 0.48 x 2T after instant of maximal 

interface height. 




fluid then diverges horizontally outwards near the bottom and moves upwards in the large 
regions between the peaks. The motion shown in figure 22 at t = 0.41 x (2T), is quite 



different from that in figure[2T] The peaks of figure 21 have collapsed into wide flat craters. 
The fluid converges inwards horizontally above the peaks, then descends into the craters 
and diverges outwards horizontally just below them. Figure 23 at t = 0.73 x (2T), shows 
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Figure 19. Temporal evolution of the amplitudes of the spatial modes with wavenumbers k c , 
2k c and v3fc c . Solid curves represent experimental results (A. Kityk & C. Wagner, private 
communication) at a w 38.5 ms -2 . Dashed curves represent the simulation for a = 38.0 ms -2 
at reso lution (in x,y,z directions) of 58 x 100 x 180. Arrows indicate times corresponding to 
figures [15}{18] 
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Figure 20. Temporal Fourier transform of the amplitudes in figure |19| Circles represent ex- 
perimental results (A. Kityk & C. Wagner, private communication) for a ~ 38.5 ms -2 . Crosses 
represent numerical results for a = 38.0 ms" 2 at resolution 58 x 100 x 180. 
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Figure 21. Velocity field at time t — 0.07 x (2T) after the instant of maximum height. Interface 
is colored according to the vertical velocity w. Arrows show velocity field at z — 0.53 mm and 
z = 6.08 mm. (Total height is 10 mm, average interface height is 1.6 mm.) For clarity, velocity 
vectors are plotted only at every fourth gridpoint in each direction. Note that the vertical and 
horizontal scales are different. One computational domain is shown. 



w (mm s" 1 ) 




Figure 22. Velocity field at time t — 0.41 x (2T) after the instant of maximum height. Vectors 
shown at z = 0.083 mm and z = 6.25 mm. Vector and color scales differ from those of figure [2l| 
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Figure 23. Velocity field at time t = 0.73 x (2T) after the instant of maximum height. Arrows 
show velocity field at z = 1.58 mm. Vector and color scales differ from those of figures [2T| and 
[221 



that the rims of the wide flat craters seen in figure [22] have in turn collapsed inwards, 
forming circular waves which invade the craters, whose remnants are visible as dimples. 
The velocity field of figure [23] shows fluid converging horizontally below these dimples. 
These are erupting at velocities which are the largest in the cycle, and will eventually 



reconstitute the high peaks seen in figure 21 



Figures [2T| [23] as well as figures [15] [18] show that these cases pose great computational 
difficulties. The interface periodically forms a very thin film (approximately 0.1 mm; see 
wide crater in figures 17 and 22) over large portions of the lower boundary, within which 
the velocity may be significant. These features make it difficult to adequately resolve the 
flow in this layer. For the time being we use a uniform grid spacing; however, in this 
case an adaptive grid would be more efficient and is under development. We have also 
simulated hexagons with resolutions of 70 x 120 x 100 and 70 x 120 x 50. Although we 
do not show these, the two highest resolutions lead to very similar spatial spectra with 
a maximum difference in amplitudes of the principal modes of about 5% between the 
70 x 120 x 100 and 58 x 100 x 180 resolutions. The case resolved by only 50 cells in the 
z direction shows differences mainly in the 2k c mode where the difference between the 
70 x 120 x 50 and 58 x 100 x 180 resolution is about 25%; for the two other modes the 
difference is about 10%. Hexagonal motifs were observed for all of the resolutions. 

The calculation for the hexagon case, for the resolution of 58 x 100 x 180 takes about 
7 h per subharmonic oscillation on a 2.16 GHz Intel processor. This corresponds to 42 h 
of calculation time for 1 s of physical time. 

In contrast to the square patterns, all of the hexagonal patterns that we have observed 
are transient. In our calculations, they last for several seconds, i.e. about 15-20 subhar- 
monic oscillation periods, over which time the amplitudes and periods of the principal 
modes remain constant. This is also the case for the experimental observations (A. Kityk 
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& C. Wagner, private communication), although the experimental lifetimes are longer. 
In our simulations, hexagonal patterns alternated with patterns with other symmetries, 
whose lifetimes were long (on the order of several seconds) but irregular. This behaviour 
suggests that the hexagonal state may belong to a heteroclinic orbit. A more extensive 
examination of the hexagonal regime will be the subject of a future investigation. 



6. Conclusion 

We have carried out full nonlinear three-dimensional simulations of Faraday waves. The 
incompressible Navier-Stokes equations for two fluid layers of different densities and vis- 
cosities are solved using a finite-difference method. The interface motion and surface ten- 
sion are treated using a front-tracking/immcrsed-boundary technique. The simulations 
are validated in several ways. First, for small oscillation amplitudes, our computations 
match the solution of Kumar & Tuckerman ( 1994[ ) to the Floquet problem which results 
from the linearized evolution equations. The boundaries of the instability tongues, i.e. 
the critical amplitude as a function of horizontal wavenumber are calculated for several 
wavenumbers on several tongues and are in good agreement with the theoretical values. 
The temporal dependence of the Floquet modes is also well reproduced by our numerical 
results, an even more quantitatively significant validation. 

For finite oscillation amplitude, our computations reproduce the square and hexag- 
onal patterns observed by Kityk et al. (2005, 2009) at moderate and high-oscillation 
amplitudes, respectively. Although the domains shown in figure [13] were chosen to ac- 
commodate square and hexagonal patterns respectively, we consider the emergence of 
these patterns at the appropriate parameter values a non-trivial test of our program, 
since these domains can also accommodate rectangles and stripes. Quantitative compar- 
isons were made between experiment and simulation of the spatio-temporal spectra. Our 
numerical results lie well within the experimental uncertainty. The hexagonal patterns 
are long-lived transients and show intriguing dynamical behavior. Our direct numerical 
simulations provide velocity fields and pressure throughout the entire domain of calcu- 
lation. Thus, we have been able to ascertain precisely the fluid motion for the Faraday 
waves, both above and below the interface between the two fluids. 

Our future studies of Faraday waves will include a more detailed investigation of the 
dynamics of the hexagonal patterns, and the simulation and interpretation of oscillons. 
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